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ABSTRACT 

In-Situ  measurements  of  the  speed  of  sound  in  the  upper  ocean  have 
revealed  the  existence  of  significant  dispersion  and  large  fluctuations 
over  the  frequency  range  25-80  kHz.  The  near-surface  values  of  c  ranged 
from  +6  M/sec  to  -3  M/sec  relative  to  the  bubble-free  value,  with  a 
maximum  estimated  error  of  0.5  M/sec.  It  was  possible  to  identify 
bubbles  of  "surface"  radius  centered  around  54  microns  down  to  4.3  meter 
depth  as  well  as  a  population  centered  around  124  microns  (at  4.3  M) 
found  at  all  depths.  The  speed  fluctuation  showed  near-Gaussian 
probability  density  functions  except  at  the  dispersion  center  frequen- 
cies. The  standard  deviation  of  the  speed  varied  from  0.27  M/sec  for 
58.0  kHz  to  0.52  M/sec  at  69.6  kHz.  The  phase  remained  temporally 
correlated  only  for  1.46  sec,  whfch  was  very  close  to  the  correlation 
times  of  the  oceanographic  variables.  The  power  spectral  density  of 
the  varying  phase  showed  its  strongest  values  at  frequencies  less  than 
0.5  Hz.;  this  was  presumably  due  to  bubbles  entrained  at  orbital  frequen- 
cies associated  with  the  surface  wave  system  which  had  maximum  energy 
in  the  same  frequency  range. 
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I.  ^INTRODUCTION 

The  objective  of  this  research  was  to  take  in-situ  acoustical 
measurements  in  the  surface  layer  of  the  ocean  in  order  to  determine 
the  existence  of  dispersion  and  fluctuations  of  the  speed  of  sound. 
The  anticipated  source  of  frequency-dependent  speed  of  sound  in  sea-water 
are  bubbles,  of  various  radii,  of  numbers  dependent  on  the  depth.  A 
knowledge  of  this  dispersive  behaviour  is  of  importance  for  sound 
propagation  in  near-surface  ducts  where  a  change  in  the  speed  of  sound 
with  depth  would  result  in  a  change  in  the  refraction  of  the  sound  ray. 
If  bubbles  exist  near  the  sea  surface  a  sound  ray  would  bounce  more 
often  from  the  surface  over  a  constant  horizontal  distance  as  compared 
with  the  bubble-free  case. 

The  problem  has  already  been  treated  in  theory  as  in  the  analysis 
worked  out  by  Meyer  and  Skudrzyk  (1950).  Fox,  Curley  and  Larson  have 
examined  in  the  laboratory  environment  the  phase  velocity  in  water 
containing  a  high  density  of  air  bubbles  0954) .  Glotov  et  al .  (1962) 
studied  a  population  of  bubbles  generated  by  wind-agitated  breaking 
waves  in  a  laboratory  tank.  Of  the  few  known  data  actually  taken  in 
the  ocean  those  given  by  Buxcey  et  al .  (1965)  and  analyzed  in  the  open 
literature  by  Medwin  (1970)  will  be  referenced  here. 

Since  virtually  no  in-situ  ocean  measurements  of  bubble  effects  were 
available,  this  research  was  proposed  and  carried  out  using  the  U.  S. 
Navy  Undersea  Research  and  Development  Center  (NUC)  Oceanographic 
Research  Tower  in  Mission  Bay/San  Diego  as  a  test  platform. 

A  stable,  accurately  known,  continuous  wave  acoustical  signal  was 
received  by  two  spatially  separated  hydrophones  mounted  on  the  sound 


axis  and  their  output  was  compared  in  phase  by  means  of  a  digital  phase 
meter.  The  phase  difference  between  the  two  hydrophones  was  recorded 
at  each  frequency.  The  frequency  was  incremented  in  steps  of  approxi- 
mately two  or  four  kHz  over  a  range  from  25  to  80  kHz  and  for  four 
different  water  depths.  The  average  results  have  been  interpreted  in 
terms  of  a  change  in  the  speed  of  sound  with  frequency  at  a  constant 
depth,  and  with  depth  at  a  constant  frequency.  The  statistics  of  the 
phase  change  with  time  have  also  been  calculated. 

At  the  same  time,  and  at  the  same  location,  the  local  sound  amplitude 
modulation  due  to  the  spatial  and  temporal  changes  in  the  medium  has 
been  examined  (thesis  by  LCDR  W.  J.  Smith,  December  1971).  In  addition, 
field  measurements  were  taken  to  determine  the  fluctuating  temperature, 
salinity,  turbulent  velocity  in  x,  y,  and  z-directions ,  speed  of  sound 
as  given  by  a  velocimeter  and  wave  height  in  order  to  allow  a  better 
judgment  of  the  acoustic-oceanic  interaction  (theses  by  Lt.  M.  Bordy, 
LCDR  C.  Duchock  and  Lt.  H.  Seymour,  March  1972), 


II.   INSTRUMENTATION 

A  steel-pipe  frame  shaped  like  a  "square  C"  of  height  and  length 
six  ft  was  used  to  mount  the  transducer  and  the  hydrophones  as  shown 
in  Fig.  1.  To  prevent  the  frame  from  acting  like  a  huge  tuning  fork 
and  thereby  displacing  the  hydrophones,  the  frame  was  put  under  tension 
halfway  down  and  at  the  open  end  by  spring-loaded  wires.  To  minimize 
the  reflections  from  the  frame  caused  by  side-lobes  of  the  transducer 
the  upper  and  lower  members  were  covered  with  acoustic  absorbent  rubber 
(SoAB)  at  the  specular  points. 

Opposite  to  the  open  end  of  the  frame  and  at  the  center  an  USRD  F27 
unidirectional  transducer  was  mounted.  At  the  sound  axis  of  this  trans- 
ducer an  Atlantic  Research  LC-10  hydrophone  was  placed  78  cm  from  the 
first  one  and  a  total  of  173  cm  from  the  F27.  This  distance  between  the 
hydrophones  was  chosen  to  make  the  measuring  path  long  enough  to  measure 
the  effect  of  the  bubbles  while  at  the  same  time  limiting  the  overall 
size  of  the  frame. 

The  transmitting  part  of  the  instrumentation  (Fig.  2)  consisted  of 
the  highly  stable  General  Radio  Coherent  Decade  Frequency  Synthesizer 
Type  1162-A  generating  a  sinusoidal  voltage  of  two  volts  rms  and  variable 
in  frequency.  This  signal  was  then  amplified  by  a  Hewlett  Packard  Power 
467A  Amplifier  to  14,1  V  rms  and  then  impressed  across  the  transducer 
emitting  the  sound  wave.  The  receiving  part  of  the  instrumentation 
(see  Fig.  2)  consisted  of  hydrophone  #1  with  a  sensitivity  of  -125  dB 
re  1  volt/microbar  and  hydrophone  #2  with  a  sensitivity  of  -106  dB  re 
1  volt/microbar  at  40  kHz  which  picked  up  the  acoustical  signal.  The 


Mounting  of  Transducer  and  Hydrophones 
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two  signals  were  amplified  by  30  dB  by  an  MUS  Pre-Amplifier  Model 
2010-030  before  being  sent  through  (150)  ft  of  waterproof  shielded 
cable.  Both  signals  were  then  bandpass-filtered  by  a  digitally  tuned 
variable  Krohn-Hite  Model  3322  Filter.  The  received  signal  was  always 
bandpass-filtered  with  the  center  frequency  equal  to  the  frequency  of 
the  transmitted  signal  through  a  passband  of  t  300  Hz,  with  -48  dB  per 
octave  attenuation  on  either  side  of  the  passband.  After  the  Bandpass- 
Filter  the  signals  were  fed  into  a  Wiltron  Digital  Phase  Meter  Model  355 
which  measures  the  phase  angle  between  the  two  ac  voltages  of  the  same 
frequency  and  provides  a  visual  display  as  well  as  a  dc  voltage  output 
proportional  to  the  measured  phase  angle  at  10  mV  per  degree.  This  dc 
voltage  was  then  FM  recorded  on  a  Sangamo  Model  3500  fourteen  track 
magnetic  tape  recorder  simultaneously  with  the  oceanographic  data  at 
a  recording  speed  of  1  7/8  inch  pe.r  sec.  The  dynamic  range  of  the 
recorder  utilized  was  .01  to  2.5  volts,  with  all  the  readings  falling 
within  this  range. 


III.  OCEANOGRAPHIC  COMPONENTS 

It  was  desired  to  define  the  volume  of  water  carrying  the  acoustical 
signal  in  terms  of  the  measurable  parameters  temperature,  salinity, 
particle  velocity,  speed  of  sound  as  given  by  a  velocimeter,  and  surface 
wave  height.  Therefore  ocean  sensors  were  mounted  close  to  the  acoustic 
system  but  not  interfering  with  its  performance  (see  Fig.  3  for  geometry) 
Three  thermistors  (Tl ,  T2,  T3)  measured  the  temporal  variation  in  the 
temperature,  and  the  temperature  probe  (T4)  of  the  Bissett-Berman 
apparatus  gave  absolute  readings  between  14°  and  19°  C.  The  Bissett- 
Berman  salinometer  yielded  the  temporal  variation  in  the  salinity  (S) 
for  values  ranging  from  37.5  to  39.5  ppt.  The  speed  of  sound  c  = 
c(  T,  S,  depth)  was  taken  by  a  Ramsey  SVTD  probe  that  provided  a  temporal 
variation  around  a  DC-value.  The  particale  velocity  was  obtained  with 
an  electromagnetic  flowmeter  reading  one  horizontal  component  (u)  or  (v) 
and  the  vertical  component  (w).  A  pressure  wave  height  indicator  and 
a  Baylor  gauge  gave  information  on  the  instantaneous  wave  height. 
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IV.  EXPERIMENTAL  SETUP 

A  continuous  wave  method  was  used  Instead  of  a  pulse-echo  system 
in  order  to  avoid  possible  errors  related  to  pulse  rise  time,  definition 
of  the  leading  edge  of  the  pulse,  distortion  of  the  pulse  when  passing 
through  the  medium  because  of  frequency-dependent  attenuation  effects 
and  the  difficulty  to  determine  the  phase  shift  on  reflection  (see  Ref. 
7),  This  concern  was  particularly  true  since  we  wanted  measurements 
over  a  wide  frequency  range. 

The  method  used  to  measure  the  speed  of  sound  is  actually  derived 
from  a  phase  interference  approach.  The  frequency  of  the  signal  applied 
to  the  transducer  is  adjusted  by  means  of  the  frequency  synthesizer 
until  a  wanted  phase  angle  is  approached  at  the  output  of  the  digital 
phase  meter.  Since  the  Sangamo  re'corder  required  a  signal  of  at  least 
100  mV  and  at  the  same  time  the  digital  phase  meter  is  most  accurate  in 
its  reading  for  phase  angles  close  to  zero,  the  frequency  was  always 
adjusted  to  get  a  phase  value  of  approximately  20  degrees.  Since  the 
digital  phase  meter  compares,  the  signal  at  input  B  with  respect  to  the 
signal  at  input  A  and  cannot  detect  integer  number  of  wave  lengths 
difference  between  the  two  signals,  the  displayed  phase  determines  the 
fractional  wave  length  by  which  signal  B  leads  or  lags  signal  A  depending 
on  the  sign.  This  means  that  an  integer  number  of  wave  lengths  plus 
this  fractional  part  just  fits  into  the  separation  distance  between  the 
two  hydrophones.  Further  increasing  the  frequency  will  finally  change 
the  phase  by  360  degrees  or  step  up  the  integer  number  of  wave  lengths 
by  one  and  thus  produce  the  next  frequency  to  be  recorded. 


Because  bubble  populations  were  expected  to  be  large  in  the  radius 
range  40-150  microns  we  examined  the  frequency  range  from  25-80  kHz. 
The  78.  cm  separation  between  the  two  hydrophones  caused  the  number 
of  wave  lengths  within  this  spacing  to  be  incremented  by  one  approxi- 
mately every   two  kHz.  In  order  to  exclude  possible  errors  in  the  results 
due  to  temporal  changes  in  the  medium  it  was  necessary  to  examine  the 
total  frequency  range  in  as  short  a  period  of  time  as  possible.  There- 
fore it  was  decided  to  accept  20  minutes  as  an  upper  limit  for 
stationarity  of  the  medium.  The  time  was  also  the  lower  limit  to  find 
the  frequency  spectrum  of  the  ocean  wave  height  and  turbulent  velocity. 
This  time  allowed  us  to  examine  the  frequency  range  in  steps  of  four 
kHz.  and  to  record  the  time-variable  phase  difference  between  the  hydro- 
phones for  one  minute  at  those  frequencies. 

In  addition  to  examining  the  speed  of  sound  with  changing  signal 
frequency  at  a  constant  depth,  readings  were  taken  of  the  variation  of 
the  speed  for  sound  with  varying  depths  at  a  constant  frequency.  Four 
different  depths  were  chosen  in  order  to  observe  near-surface  to  near- 
bottom  effects:  4.3,  9.3,  14.3,  7.3  m.  These  depths  were  chosen  to 
minimize  interference  from  the  underwater  beam  structure  of  the  tower. 
The  order  was  picked  to  decrease  and  to  check  on  the  possibility  of 
longtime  variations  in  the  medium.  The  depths  are  in  meters  below  mean 
water,  where  mean  water  was  6.7  m  below  the  concrete  deck  of  the  NUC  tower. 

All  the  readings  were  taken  within  4  1/2  consecutive  hours  at  night 
and  early  in  the  morning  on  October  22,  1971.  The  obtained  data  are 
given  in  Appendix  A. 


V.  DATA  REDUCTION 

Five  twenty-minute  runs  (split  into  ten  two-minute  individual 
records)  make  up  the  recorded  analog  data.  The  approach  used  was  to 
read  this  information  out  on  a  strip-chart  recorder,  to  digitize  the 
printed  analog  data,  and  to  use  a  computer  to  carry  out  a  statistical 
and  spectral  analysis. 

The  phase  fluctuation  and  the  other  oceanographic  parameters  were 
initially  recorded  at  1  7/8  inches  per  second.  The  tape  recorded  data 
were  played  back  at  7  1/2  inches  per  second  and  first  recorded  on  an 
eight  channel  Brush  Mark  200  strip-chart  recorder  in  order  to  allow 
visual  comparison  between  the  various  sensor  outputs  (Fig.  4). 

The  tape  recorded  data  was  then  played  back  a  second  time  at  7  1/2 
inches  and  recorded  on  a  two  channel  Brush  strip-chart  recorder.  This 
compressed  the  20  minute  run  to  five  minutes  and  resulted  in  an  increase 
in  frequency  by  a  factor  of  four.  The  strip-chart  recorder  was  set  at 
25  mm/sec.  chart  speed  and  10  mV/div  sensitivity,  corresponding  to 
1  degree/div. 

The  strip-chart  records  for  each  frequency  at  each  depth  were  digitized 
using  the  Fleet  Numerical  Weather  Central  facility  tracing  digitizer  at 
Point  Pinos,  California.  The  digitizer  sampled  the  trace  in  increments 
of  hundredths  of  an  inch  in  the  y-direction  for  each  one  one-hundredth 
of  an  inch  advance  in  the  x-direction  and  recorded  this  on  a  nine-track 
magnetic  tape.  Every  record  was  digitized  by  first  following  the  zero 
Volt  reference  line  on  the  strip-chart  and  then  moving  the  tracer 
vertically  at  the  start  of  a  record  from  the  base-line  onto  the  trace. 
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By  doing  this  it  was  assured  that  the  digitized  data  points  were 
absolute  values  referenced  to  zero  Volt.  Thus,  the  scaling  factor  in 
the  y-direction  is 

.254  mm/data  point  -  ,254  degree/data  point 
resulting  in 

.254  mm/d,p.  X  jr-  sec/mm  =  .01016  brush  sec/d.p. 

=  ,04064  real  sec. /data  point 

This  implies  that  the  sampling  rate  is  25  Hz  and  the  Nyquist  rate  is 
then  12.5  Hz. 

The  digitized  data  were  read  out  from  the  magnetic  tape  and  punched 
on  data  cards  using  the  Fleet  Numerical  Weather  Central's  CDC  5000 
computer  (see  Computer  Program  1).  The  data  cards  could  then  be  used 
for  analysis  on  the  USNPGS  IBM  360  digital  computer. 


VI.  COMPUTATIONAL  METHOD 

Since  we  were  working  with  a  continuous  wave  method  the  signal 
received  by  the  first  hydrophone  can  be  written  as 

yA  =  a  sin  (cot) 

and  the  signal  received  by  the  second  hydrophone  B  has  then  to  be  written 

y„  =  b  sin  (wt-kx) 

where     k  =  the  wave  number  =  -^—  =  —^— 

c      c 

f  =  the  frequency  of  the  applied  signal 

c  =  the  phase  wave  velocity 

x  =  the  separation  between  the  hydrophones 

Since  we  compare  the  two  signals  with  respect  to  phase  by  using  the 
digital  phase  meter  we  have  to  find  values  for  w  for  which 

kx  =  \p 
where  ij;  is  the  total  phase  difference  between  the  two  signals  due  to 
an  integral  number  of  wave  lengths  and  a  fraction  of  the  wave  length. 
We  write 

kx  =  (2ir}n  +  <(>        n  =  an  integer 

kx    =     ~x     =     (2^)n  +  <f> 


0)X 

c  C2iT)n  +  <)> 


fx 

1  -    ' •  where  $  is  in  radians 


n  +  W11- 

£7T 


The  final  result  will  be 


fx 
c  =  where  $   is  in  degrees 

n  +  ^ 


360 


Therefore  the  speed  of  sound  may  be  obtained  by  using  the  measured 
frequency  at  which  a  phase  difference  <$>   was  noted.  In  addition,  x  and 
n  must  be  known.  The  integer  n  can  be  found  by  calculating  c  from  the 
Wilson  equation  (Kinsler  and  Frey,  Eq.  15.1)  using  readings  for  the 
temperature  and  the  salinity.  We  find 


% 

a 

C 

f 

and  then 

n 

= 

X 

'A 

VII.  BUBBLE  THEORY 
The  speed  of  sound  through  a  fluid  is  given  by 

c  =  ~T~ 

where    w  =  sound  frequency 

k  =  2tt/a  =  propagation  constant  in  bubble-free  water 

Following  the  analysis  carried  out  by  Morse  and  Feshbach,  an  incident 
acoustic  pressure,  of  wave  length  yery   large  compared  to  the  bubble 
radius  R,  having  q  velocity  potential  4>.e~  w  ,  will  produce  a  total 
field  averaged  over  all  possible  sizes  and  configurations  of  bubbles 
when  written  in  differential  equation  formulation  of 

v2<  O  +  kh2(r)  <*>  =  0 


where 

khCr)     *     k    +    &  /    nirALlJR 

D         K  \    (u)o/u>r  -  l  -  i$0ooA>) 

r  =  coordinate  position  of  scattering  region 

k.  =  propagation  constant  in  bubbly  water 

R  ■  bubble  radius 

w  =  resonant  frequency  of  bubble 

$  =  damping  constant  of  bubble 

n(R)dR  =  number  of  bubbles  per  unit  volume  in  radius  increment 
R  to  R  +  dR 

The  quantity    l  l'  '   1S  the  index  of  refraction  in  the  bubbly  region 


ond  the  real  part  defines  the  local  speed  of  sound.  This  term  is 

,2 


Re(kh(r)}  -  k  +  f!  /   B  nffi)  [M*  '  U   <"* 


r«y 


c  =  k  = 


+  ^fR  n^R)  i>0/a))2 : 1]  dR 

k/   .[(co    /co)2    -   l]2   +6  2(0)    /co  )2 
R  °  ° 


L     2,c0  <:  Ri  "(y^)-i]dRt 

c  ~*~2~~    ^ 

1      LCwQ/aj)      -  1  j     +  6    (<oq/co) 


where 


nCR)  dR    =    iiiEMR 

4/3  ttRj 

aCR,)  dR. 

n(R  )dR    =    } _L 

1        1         4/3  uR.J 

uCR)  dR    ■    fractional   air  to  water  volume  for  the  increment  between 
R  and  R  +  dR 

u(R.)dR.  =     fractional   air  to  water  volume  for  the  increment  between 
1       1         R.  and  R.  +  1/n 

Therefore  we  can  write 

1 


c    = 


JL,    +     3cQ  uCRt)  dR.         [(a)0j/co)2  :  1] 

co        7T  C~     n 


2w  .  Ri  [(co     /co)2   ~   l]2   +    S.2(co     /co)2 

1  LV   Oj'     '  J  1     v    0-j       ' 

This  is  the  dispersion  relation  that  predicts  a  decrease  of  sound  speed 
caused  by  bubbles  of  resonance  frequency  greater  than  the  incident 
frequency  and  an  increase  of  sound  speed  for  bubbles  with  lower  resonant 
frequencies  than  the  impinging  sound. 


The  resonance  frequency  of  clean  air  bubbles  in  water  is  given  by 
(Eckart,  1945) 

fo  =  o)o/2tt  =  Cl/2TrRo)IC3YPo/po)Cg/a)]+1/2 

where    y  =  ratio  of  specific  heats  of  bubble  gas  (=1.402  for  air) 
P  -  ambient  pressure  at  bubble  depth 
p  =  density  of  surrounding  water 

The  factor  g  is  due  to  the  influence  of  surface  tension,  it  differs  from 

unity  only  for  s/ery   small  bubbles,  for  the  smallest  bubbles  of  this 

study  (R  =  50  microns)  g  =  1.02.  The  factor  a   modifies  the  value 

of  y,  i.e.,  ycT  is  the  effective  ratio  of  specific  heats  and  lies  for 

the  smallest  bubbles  of  this  study  in  the  range  1,24  <  y«~  ^_1.40. 

The  two  effects  are  to  some  extent  mutually  counteracting  in  their 

effect  on  the  resonance  frequency  and  therefore  the  term  yg/a  will 

be  assumed  to  equal  1.33. 

The  damping  factor  as  a  function  of  bubble  resonance  frequency 

as  given  by  Devin  will  be  approximated  by 

0,322 
b\     -     <027Hk/2,  x!03^ 


VIII.  EXPERIMENTAL  RESULTS 

Five  different  freqency  sweeps  have  been  carried  out,  runs  6  and  7 
at  a  depth  of  4.3  m  and  the  remaining  three,  runs  8,  9,  and  10,  at 
9.3,  14.3  and  7,3  m  depth,  respectively. 

Run  6  covered  a  total  frequency  range  from  25  ~  80  kHz  in  steps  of 
approximately  two  kHz.  to  be  able  to  detect  the  detailed  variations 
in  the  speed  of  sound;  for  the  remaining  runs  readings  were  taken  every 
four  kHz.  in  order  to  examine  the  frequencies  of  interest  in  as  short 
a  period  of  time  as  possible.  Run  7  started  on  hour  after  the  end  of 
run  5,  but  for  the  five  runs  reported  here  all  the  readings  were  taken 
within  4  1/2  consecutive  hours  between  0350  and  0820  on  October  22,  1971. 

Single  temperature  and  salinity  readings  taken  at  each  depth  and 
a  single  bathythermograph  record  taken  at  8  a.m.  (Fig.  5)  allow  one 
to  calculate  the  speed  of  sound  for  bubble-free  water;  the  Wilson 
equation  is  used  in  the  following  form: 

c  =  1449  +  4,6t  -  0.055t2  +0.0003t3  +  (1  ,39-*,0l2t)  (S-35)  +  0.017d 

where    c  =  speed  of  sound  in  meters  per  second 

t  =  temperature  in  degrees  centigrade 

S  =  salinity  in  parts  per  thousand 

d  =  depth  in  meters 
The  bubble-free  values  obtained  for  4.3,  9.3  and  14.3  m  are  used  as 
references  for  the  dispersion  curves  and  are   included  in  Fig.  6. 

The  calculated  values  for  the  speed  of  sound  based  on  the  obtained 
data  are  examined  and  analyzed  with  respect  to  a  frequency  —  and/or  depth 
dependence.  If  there  is  any  change  in  c  with  frequency  then  this  change 
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Fig.  5. 

Bathythermograph  Approximation 
for  22  October  1971  ,  0800 
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is  attributed  in  general  to  either  molecular  relaxation  or  bubbles 
since  for  the  sea-water  medium  no  other  frequency-dependent  mechanisms 
are  known  at  the  frequencies  of  this  study,  But  the  magnitude  of  the 
molecular  relaxation  effect  can  be  shown  to  be  negligible  (a  frequency 
shift  from  20  to  70  kHz.  results  in  a  Ac ^0.3  m/sec.)  as  compared  to 
the  bubble  effects. 

A,   FREQUENCY  BEHAVIOUR 

All  the  calculated  values  of  c  as  a  function  of  frequency  for  the 
different  runs  have  been  plotted  on  Fig.  6  for  comparison.  They  should 
be  interpreted  by  making  use  of  the  reference  value  for  the  non-bubbly 
medium  which,  because  of  the  temperature  change,  shifts  by  as  much  as 
4.2  m/sec.  in  going  from  4.3  m  to  a  depth  of  14.3  m. 

The  first  point  of  interest  is  the  dispersive  variation  of  c  at  a 
single  depth.  The  first  run  (6)  was  carried  out  at  a  depth  of  4.3  m 
and  examined  frequencies  from  30.7  up  to  79.2  kHz  during  a  40  minute 
period  of  time.  Most  significant  is  the  strong  increase  in  the  speed 
of  sound  by  as  much  as  6  m/sec.  compared  to  the  non-bubbly  value  for 
c  of  1510.47  m/sec,  and  caused  by  the  peak  at  34.7  kHz  followed  by  a 
smaller  negative  change  of  3  m/sec  at  40.6  kHz.  For  frequencies  from 
47  to  60  kHz  the  speed  of  sound  shows  nearly  no  dispersive  behaviour. 
The  higher  frequencies  are  characterized  by  another  peak  in  the  c-curve 
roughly  one  fourth  the  strength  of  the  peak  at  34.7  kHz. 

Run  seven  is  interesting  since  it  covered  the  same  frequency  range 
and  the  same  depth  as  run  six  but  for  the  major  part  now  in  steps  of 
four  kHz  rather  than  2  kHz.  Since  this  allows  a  comparison  between  the 
two  runs,  a  reminder  is  necessary  that  there  is  a  time  lag  of  one  hour 
between  the  end  of  run  six  and  the  start  of  run  seven.  This  time  delay 


may  account  for  the  slight  shift  in  the  dispersion  curve  downward  and 
to  the  right.  The  shape  matches  the  previous  one  pretty  well  ,  with  two 
peaks  clearly  discernible. 

In  run  eight  the  depth  is  increased  by  five  meters  as  compared  with 
the  previous  run.  The  main  feature  is  a  resonance  pattern  with  a  weaker 
maximum  at  36.8  kHz  rather  than  34.7  kHz  at  shallower  depth. 

Run  nine,  at  the  greatest  depth  (J4.3  m) ,  differs  in  the  pattern 
for  the  speed  of  sound  over  the  frequency  range  by   showing  an  increase 
in  c  over  a  wide  frequency  band  with  a  maximum  at  approximately  41  kHz 
and  no  second  peak.  The  obvious  shift  in  the  level  is  caused  by  a 
decrease  of  the  non-dispersive  value  of  c  from  1510.47  at  4.3  m  to 
1506.27  m/sec  at  14.3  m. 

Run  ten  was  made  back  at  a  depth  of  7.3  m.  The  two  readings  at 
57,7  and  61.59  kHz  are  believed  to  be  non-representative  since  a  boat 
approaching  the  tower  injected  a  large  number  of  bubbles  into  the  volume 
under  consideration.  Excluding  these  data  points,  the  dispersion  curve 
shows  a  main  peak  around  35.5  kHz. 

Summarizing  the  dispersion  data  it  can  be  stated  that: 

1)  The  change  of  the  speed  of  sound  with  frequency  at  sea  was  observed 
to  be  significantly  stronger  than  expected. 

2)  The  change  in  the  speed  of  sound  does  not  follow  the  usual  dispersive 
pattern  for  a  single  dominant  bubble  regime. 

3)  The  speed  of  sound  curve  at  4,3  m  has  two  peaks.  The  frequency  of 
the  lower  peak  changes  from  34.7  kHz  at  4.3  m  to  41  kHz  at  14.3  m,  which 
is  a  shift  toward  higher  frequencies  as  expected  for  surface-generated 
bubbles  carried  downwards  "isothermally"  by  turbulence.  Hov/ever,  the 
shift  is  smaller  than  expected:  isothermally  convected  identifiable, 


non-diffusing,  surface-initiated  bubbles  which  resonate  at  35  kHz  at 
4.3  m  would  resonate  at  54  kHz  and  not  at  the  observed  41  kHz  when 
they  reach  14,3  m.  There  are  two  explanations  that  can  be  proposed 
for  this  discrepancy: 

a)  the  presence  of  the  higher  frequency  peak  in  the  near-surface 
measurement  has  shifted  the  low  frequency  peak  to  higher  frequencies 
than  it  would  have  had  by  itself.  The  absence  of  the  second  peak  in 
the  deeper  measurement  (due  to  more  rapid  disappearance  of  smaller 
bubbles)  permits  the  low  frequency  peak  to  more  accurately  identify 
the  resonance  frequency  of  these  bubbles.  We  herein  assume  that  near 
the  surface  we  encountered  a  bubble  population  made  up  of  larger  surface- 
generated  bubbles  and  smaller  bubbles  introduced  primarily  by  continental 
aerosols  that  drop  into  the  water.  To  check  on  the  dispersion  curve  that 
results  from  a  specific  bubble  distribution,  a  bubble  population  was 
created  peaking  around  60  and  120  microns  as  given  in  Fig.  7.  Subse- 
quently the  computer  program  (3)  was  written  calculating  the  speed  of 
sound  dispersion  based  on  the  equations  derived  in  Section  VII  and  using 
the  above  mentioned  volume  concentration  of  bubbles.  The  result  is  a 
dispersion  curve  given  in  Fig.  8  which  is  in  good  agreement  with  our 
data  obtained  at  4.3  m.  The  observed  dispersion  curves  for  the  other 
depths  support  this  hypothesis,  since  only  part  of  the  larger  bubbles 
are  carried  down,  decreasing  in  number  with  increasing  depth  whereas  the 
smaller  bubbles  rapidly  disappear  due  to  gas  diffusion  out  of  the 
bubbles. 

b)  the  larger  bubbles  originate  in  the  volume  or  at  the  bottom,  rather 
than  at  the  surface.  As  they  rise  they  lose  gas  by  diffusion,  so  that 
when  they  reach  the  surface  they  are  not  as  large  as  predicted  by  the 
non-diffusion  assumption  and  are  thereby  at  higher  frequencies  than  expected. 
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B.  DEPTH  BEHAVIOUR 

Next  the  variation  in  c  with  depth  for  a  constant  frequency  will 
be  considered,  On  the  assumption  that  the  bubble  densities  were  not 
a  function  of  time  during  the  2  1/2  hours  of  runs  7-10,  all  these  data 
have  been  rearranged  to  display  the  variation  in  the  speed  of  sound 
with  depth  at  various  constant  frequencies  as  given  in  Table  I,  These 
values  have  been  used  to  plot  the  speed  of  sound  versus  depth  in  Fig,  9, 
These  straight  line  segments  reflect  the  major  variations  in  c  with 
frequency  as  described  in  the  last  section. 

In  order  to  examine  the  variation  of  the  speed  of  sound  versus 
depth,  it  is  necessary  to  subtract  out  first  all  possible  temperature, 
salinity,  and  pressure-based  changes  with  depth.  To  do  this,  use  is 
made  of  the  following  coefficients  for  variation  with 

temperature  jj    =  +  5  ft/ (sec)  (°F) 

salinity  ||  =  +  4  ft/ (sec)  (ppt) 

depth  ||  =  +  0.017  ft/Csec)  (ft) 

The  salinity  changed  by  -.1  ppt  over  10  m  and  therefore  accounts  for  a 
change  of  c  with  S  of  -.4  meter  per  second.  It  is  necessary  to  add  to 
this  a  change  of  +.17  m/sec  due  to  the  depth  effect  such  that  the  total 
effect  of  salinity  and  depth  on  c  will  be  -.23  m/sec  when  going  from 
4,3  to  14.3  m;  this  small  effect  will  be  neglected. 

The  temperature-caused  changes  in  the  speed  of  sound  have  to  be 
taken  into  consideration.  The  BT  curve  sketched  in  Fig.  5  is  approxi- 
mated by  straight  line  elements  and  is  used  with  the  Wilson  equation 
to  calculate  values  for  the  speed  of  sound  at  different  depths  for  the 


TABLE  I.  Speed  of  Sound  in  m/sec 


^Depth  (m) 

4.3 

4.3 

7.3 

9.3 

14.3 

Freq.    (kHz) 

28.97 

1505.73 

1503.28 

30.70 

1511.99 

1511.37 

1507.67 

32.77 

1514.18 

1513.88 

1509.89 

34.76 

1516.46 

1516.94 

1512.51 

36.79 

1512.90 

1516.52 

1511.55 

1509.2 

38.73 

1509.48 

1511.82 

1510.93 

40.61 

1507.30 

1510.22 

1509.89 

42.50 

1508.13 

1509.30 

1509.53 

44.47 

1509.39 

1510.17 

1509.89 

46,40 

1510.29 

1509.62 

48.36 

1511.31 

1510.05 

1509.84 

50.27 

1511.36 

1508.78 

52,20 

1511,03 

1509.79 

1509.36 

54,12 

1510.74 

1509,20 

56.03 

1510.92 

1509.99 

1510.66 

1506.56 

59,88 

1510.31 

1509.19 

1510.64 

1506.78 

63,67 

1509.11 

1507.65 

1508.95 

1505.43 

65.62 

1508.11 

1507.86 

67.54 

1509.14 

1507,55 

1507.92 

71.42 

1501,32 

1509.26 

1508.42 

1505.01 
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the  non-bubble  case.  After  doing  this  it  is  possible  to  compare  the 
measured  values  of  c  with  the  values  for  bubble  free  sea  water.  The 
plot  of  the  temperature-caused  changes  has  to  be  referenced  to  a  value 
of  c  that  holds  for  the  bubble-free  medium  and  the  previous  calculated 
value  of  1506,27  m/sec  at  14.3  m  was  used.  The  non-bubble  theoretical 
value  of  c  is  subtracted  from  the  measured  value  and  the  difference  is 
then  called  ac^,  due  to  bubbles,  which  is  plotted  on  the  same  graph. 
Since  the  shift  in  the  frequency  peaks  causes  the  pattern  of  the  measured 
c  versus  depth  plot  to  be  markedly  different  for  the  various  frequencies, 
five  different  frequencies  have  been  picked  that  are  felt  to  be  repre- 
sentative for  the  total  frequency  region:  36.79,  44.47,  56.03,  63.67, 
71.42  kHz. 

This  representation  of  c  versus  depth  at  a  constant  frequency  can 
only  be  used  to  examine  whether  one  encounters,  at  a  specific  depth, 
resonant  bubbles  such  that  for  sound  frequencies  smaller  than  the  bubble 
resonance  frequency  one  experiences  a  decrease  of  sound  speed  and  for 
impinging  sound  of  higher  frequencies  than  the  resonance  frequency  an 
increase  in  c.  It  should  be  clear  that  this  plot  ignores  depth  effects 
on  an  identifiable  surface-generated  bubble  carried  to  different  depths. 

Figure  10  presents  this  information  for  a  frequency  of  36.79  kHz. 
Since  there  is  a  strong  positive  acb  of  5.25  m/sec  at  4.3  m  this  means 
that  we  encounter  strong  populations  of  bubbles  with  lower  resonance 
frequencies  than  the  impinging  sound,  in  other  words,  there  are  many 
bubbles  that  resonate  below  36.79  kHz.  Figure  11  represents  the  situ- 
ation at  44.47  kHz  and  shows  that  we  sense  bubbles  of  resonance  frequency 
greater  than  the  incident  frequency  at  4.3  m  depth  and  of  resonance 
frequency  smaller  than  the  incident  frequency  at  14,3  m.  The  picture 
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changes  v/hen  going  to  Fig.  12.  At  a  frequency  of  56.03  kHz  there  is 
no  significant  dispersive  behaviour  at  all.  For  frequency  63.67  kHz, 
Fig.  13,  at  all  depths  we  sense  bubbles  of  resonance  frequencies  greater 
than  the  frequency  of  the  acoustical  signal.  Since  the  difference 
between  the  measured  value  and  the  non-bubble  theoretical  value  is 
much  larger  at  4.3  than  at  14.3  m  we  are  either  very   close  in  frequency 
to  resonant  bubbles  at  4,3  m  or  the  population  is  considerably  stronger 
than  at  14.3  m.  Figure  14,  71.43  kHz,  finally  shows  a  smaller  negative 
acb  with  only  minor  differences  between  the  various  depths. 

Summarizing,  we  can  make  the  following  statements  related  to  the 
different  depths: 

4.3  m:  There  is  a  large  population  of  bubbles  having  a  resonance 
frequency  close  to  and  somewhat  smaller  than  36.79  kHz.  At  44.47  and  at 
56.03  kHz  we  sense  and  at  63.67  kHz  we  detect  another  population  of 
bubbles  with  an  w  greater  than  63.67  kHz. 

7.3  m:  Bubbles  with  a  resonance  frequency  less  than  36.79  kHz  are  also 
detected  at  this  depth.  At  44.47  kHz  we  are  far  enough  away  from  this 
region  and  encounter  no  dispersive  characteristic. 
9.3  m:  There  is  a  population  of  bubbles  resonating  below  36.79  kHz. 
14.3  m:  There  are  bubbles  of  resonance  frequencies  below  44.47  kHz. 
At  56.03  kHz  this  region  of  resonance  frequencies  has  been  passed  and 
no  dispersive  behaviour  is  shown. 
Conclusions: 

At  4,3,  7.3,  and  9.3  m  we  observed  bubbles  of  resonance  frequency  below 
36.79  kHz  and  above  71.42  kHz  and  at  14.3  m  we  found  a  broad  region  of 
bubbles  resonating  below  44.47  kHz.  The  sound  speed  dispersion  is 
strongest  at  4,3  m  and  is  reduced  in  amplitude  when  going  to  a  greater 


depth.  The  conclusions  derived  from  the  variation  in  the  speed  of  sound 
with  depth  are  consistent  with  the  observations  made  for  the  variation 
in  the  speed  of  sound  with  frequency. 


IX,  STATISTICAL  AND  SPECTRAL  ANALYSIS 

The  signal  received  at  hydrophone  one  does  not  have  a  constant  phase 
difference  with  respect  to  that  at  hydrophone  two  but  is  a  random  signal. 
This  is  the  equivalent  of  a  noisy  signal  where  the  noise  in  this  case  is 
due  to  the  integrated  effects  of  time-varying  speed  of  sound  due  to 
varying  temperature,  salinity  and  bubble  patches  along  the  sound  path. 
It  is  assumed  that  the  random  signal  exhibits  statistical  stationarity 
which  is  partially  described  in  terms  of  a  probability  density  function 
(P.D.F.)  and  the  statistical  moments.   In  addition  to  a  statistical 
description  it  is  valuable  to  use  the  time  series  to  determine  the 
temporal  correlation  function  and  the  power  spectral  densities  of  the 
time-varying  phase  (Ref.  '2). 

The  variable  phase  signal  that  we  received  at  a  specific  frequency 
and  for  a  period  of  approximately  90  sec  is  a  continuously  changing 
random  variable  within  some  finite  interval.  But  since  this  recorded 
information  is  digitized,  the  random  variable  can  assume  only  a  finite 
number  of  distinct  values  so  that  we  must  treat  it  as  a  discrete  random 
variable.  In  practice  the  recording  technique  used  resulted  in  magni- 
tude increments  of  .254°  spaced  in  time  by  .04064  sec.  For  one  frequency 
record  we  then  get  n  data  points.  If  A  is  one  of  the  possible  data  points, 
then  n.  is  the  number  of  times  A  occurred  and  the  ratio  n./n  is  the 
relative  frequency  of  occurrence  of  the  event  A  for  that  particular 
record.  This  approach  allowed  the  use  of  the  sampled  data  points  to 
form  a  histogram  that  can  be  interpreted  as  a  distribution  function: 

FCx)   =   P(*  •:  X) 


from  which  a  density  function  can  be  obtained: 

p(x)  =  Z    P(x.)  <5(x  -  x.) 
i 

Examining  these  histograms0  it  turns  out  that  the  main  pattern  of  the 
envelope  follows  a  Gaussian  distribution.  This  can  be  expected  since 
in  the  experimental  measurement  many  irregular  and  fluctuating  causes 
like  weak  temperature  and  salinity  inhomogeneities  caused  the  phase 
reading  to  have  small  random  variations  which  produce  the  Gaussian 
P.D.F.  about  the  undisturbed  value.   It  seems  that  the  bell-shaped 
curve  experiences  stronger  changes  around  35  and  65  kHz.  This  different 
behaviour  is  pointed  out  by  a  series  of  historgrams  (Figs.  15-19)  for 
the  latter  frequency  range  and  for  4.3  m  depth.  At  59.9  kHz  the  general 
shape  is  narrow  Gaussian  but  at 63.7  kHz  the  distribution  is  nearly 
evenly  spread  over  a  wider  center  range  that  sharply  drops  at  the 
skirts.  At  67.6  kHz  the  envelope  changes  again  and  achieves  at  69.6 
kHz  a  wide  Gaussian  pattern.  At  77.3  kHz  the  pattern  is  again  that  of 
a  narrow  Gaussian  noise.  Since  this  change  in  the  density  function 
is  a  particular  phenomenon  for  frequencies  around  35  and  65  kHz  and  is 
not  found  for  frequencies  in  between,  it  is  postulated  that  time  varying 
resonant  bubble  populations  predominantly  around  frequencies  35  and  65 
kHz  caused  the  variation  with  frequency.  These  frequencies  are  essen- 
tially the  same  as  observed  previously  in  the  peaks  in  the  analysis  of 
speed  as  a  function  of  frequency. 

The  statistical  average,  or  mean  value,  x*  of  the  random  variable  <f> 
has  been  used  to  calculate  the  mean  speed  of  sound  which  is  the  expected 
value  for  c.  By  calculating  the  variance 
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*x   =  (x  -  x)   =  x  -  x 
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and  taking  the  square  root  of  the  variance  to  obtain  the  standard 
deviation,  a     ,  one  has  a  measure  of  the  spread  of  the  observed  values. 
A  large  standard  deviation  implies  a  large  spread  of  likely  values  for 
any  given  observation,  and  vice  versa.  The  standard  deviation  was 
calculated  for  all  frequency  readings  from  30.7  to  79.2  kHz  and  is 
plotted  as  a  function  of  frequency  in  Fig.  20. 

Since  there  was  a  time  lag  of  two  minutes  between  frequency  runs, 
o  is  for  this  case  not  only  a  function  of  frequency  but  also  of  time, 
which  is  increasing  with  frequency.  To  be  able  to  judge  which  fraction 
of  the  variation  in  a  is  a  time  effect  and  which  is  the  frequency 
effect  a  continuous  20  min  record  of  phase  modulation  at  60  kHz  was 
split  into  ten  two-minute  records  so  that  the  standard  deviation  could 
be  examined  as  a  function  of  time  alone  (Fig.  21).  Comparing  the  two 
graphs  there  is  a  mean  value  of  1.28°  and  a  total  variation  of  i  .32° 
for  the  time  record;  there  is  a  much  stronger  excursion  for  the  frequency 
plot  which  shows  a  mean  value  of  1.92°  at  60  kHz  with  a  deviation  of 
+  1.92°  at  40.6  and  of  +  2.56°  at  69.6  kHz.  One  can  rule  out  temperature 
and  salinity  patches  as  the  causing  agency  for  the  effect  in  Fig.  20 
since  those  would  not  show  a  frequency  dependent  spread  of  the  measured 
values.  It  seems  that  the  presence  of  bubbles  in  the  medium  causes  this 
strong  change  in  the  standard  deviation  as  a  function  of  frequency.  The 
frequencies  of  the  peak  values  of  a     are  approximately  the  same  as  the 

A 

frequencies  already  identified  as  centers  of  resonant  bubbles. 

Next  the  autocorrelation  function  RCt)  will  be  considered  since  it 
provides  information  about  4>(t)  and  also  relates  to  the  frequency-domain 
description  of  the  random  signal.  R(j)  is  defined  by 

RCx)  =  <  x(t)x(t+T)>/<x(t)x(t+T)> 
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Fig,  21.  Standard  Deviation  as  a  Function  of 
Time  Only  at  Frequency  60  kHz 


It  should  be  noted  that  the  autocorrelation  function  is  a  measure  of 
both  temporal  variation  and  statistical  dependence.  The  plots  of  the 
normalized  autocorrelation  at  30.7  and  46.4  kHz  are  considered  to  be 
representative  and  are  given  in  Figs.  22  and  23.  The  graph  for  30.7 
kHz  shows  a  variation  superimposed  on  the  almost  ideal  damped  cosine 
pattern'  obtained  at  46.4  kHz.  This  variation  is  probably  due  to  bubble 
effects.  Recall  that  the  frequency  30.7  kHz  was  observed  as  a  predo- 
minant resonance  frequency  in  the  speed  of  sound  discussion  of  section 
VIII.  1  he  time  it  takes  for  the  autocorrelation  to  drop  to  1/e  of  its 
value  is  plotted  versus  frequency  on  Fig,  24,  No  specific  regularity 
can  be  found.  On  the  average  this  "correlation  time"  is  approx.  1.8  sac. 

Examining  the  autocorrelation  function  we  arrive  at  the  following 
results:  The  phase  modulation  remains  correlated  (1/e)  for  1.46  sec, 
the  salinity  for  3.96  sec  and  the  speed  of  sound  as  measured  by  the 
velocimeter  for  2.28  sec  (Figs.  25-27).  The  temperature  records  have 
not  yet  been  analyzed  (thesis  of  LCDR  Duchock  to  be  issued  March  1972). 
Since  for  the  velocimeter  (operating  at  3  MHz)  c  =  c  (T,  S)  it  appears 
that  the  temporal  variation  in  c  is  due  principally  to  salinity  patches 
(and  most  probably  temperature  patches).  On  the  other  hand,  our  phase 
measurement  is  a  function  of  not  only  S  and  T  but  also  of  the  bubble 
distribution.  The  still  smaller  correlation  time  of  1.46  sec  is  evidence 
of  additional  decorrelation  due  to  bubbles. 

Turning  to  the  frequency  domain,  the  power  spectral  density  of  the 
phase  variation  gives  the  distribution  of  average  power  of  the  fluctu- 
ations in  frequency.  To  find  the  PSD  the  Mean  Lag  Product  Method  by 
Blackman/Tukey  was  used.  The  main  feature  of  this  method  is  to  calculate 
not  only  the  autocovariance  function  but  tn  addition  an  apparent 
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Ftg.  22.     Autocorrelation  at  30.7  kHz 
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Fig.  25.  Autocorrelation  for  Phase  Modulation  at  fj(j  kHz 
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Autocorrelation  for  Speed  of  Sound  Fluctuations 
by  Velocimeter 
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Fig.  27.  Autocorrelation  for  Salinity  Fluctuations 


40 


autocovariance  by  using  a  Parzen-lag-window.  By  taking  the  Fourier 
transform  of  this  function  one  gets  a  smoothed  power  spectrum.  The 
result  of  this  computation  using  computer  program  (2)  is  a  plot  of  phase 
angle  squared  per  Hz  versus  frequency  (The  frequency  resolution  is  .05 
Hz.  For  better  resolution  in  the  amplitude,  10  log1Q  of  the  ordinate 
is  taken  such  that  the  dimension  is  in  dB),  Several  conclusions  can 
be  derived  from  the  PSD-plot  for  the  4.3  m  depth.  For  nearly  all 
frequencies  there  is  a  distinct  peak  around  .48  Hz  probably  due  to  the 
orbital  effects  caused  by  the  observed  predominant  ocean  wave  component. 
From  0.5  to  2.5  Hz  the  values  for  the  PSD  follow  an  exponentially 
decaying  curve.  This  pattern  is  interrupted  by  peaks  mainly  at  .8, 
1,3,  1,55,  1.95,  and  2.35  Hz  that  may  show  up  at  two  or  at  most  three 
consecutive  readings  (see  peak  at  2.3  Hz  in  Figs,  28-31),  or  occur  at 
one  frequency  and  the  next  time  two  records  later.  Since  every  recording 
took  two  minutes  this  then  suggests  that  the  agency  forcing  these  changes 
in  the  frequency  spectrum  is  of  a  nature  that  can  result  in  variations 
that  show  up  only  over  a  time  period  of  four  minutes  or  repeat  ewery 
four  minutes. 

In  general  temperature,  salinity,  and/or  bubble  patches  carried  by 
turbulence  and  crossing  the  acoustical  path  could  be  the  driving  force 
for  these  variations.  Since  the  sound  signal  integrates  all  these 
effects  it  is  necessary  to  examine  the  temperature  and  salinity  records. 
For  the  20  min  run  and  the  depth  of  4.3  m  the  power  spectral  density 
of  the  phase  modulation  with  a  frequency  resolution  of  .0045  Hz,  of  the 
Ramsay  velocimeter  reading,  and  of  the  Bissett-Berman  salinometer,  all 
follow  approximately  the  same  exponentially  decaying  pattern,  Figs. 
31-33  (these  graphs  courtesy  of  LCDR  Ouchock  and  LT  Seymour).  The 
principal  difference  is  that  the  PSD  of  the  phase  shows  more  power 
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Fig.  28,  Power  Spectral  Density,  46.4  kHz 
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Fig.  29.  Power  Spectral  Density,  48.4  kHz 
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Fig.  30.  Power  Spectral  Density,  50.3  kHz 
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2.5 


concentrated  below  0.5  Hz  where  orbital  velocities  are  strong.  All 
three  plots  of  PSD  show  good  agreement  even  in  the  fine  structure  in 
the  range  from  0.5  -  1.0  Hz. 


X.  ERROR  ANALYSIS 

The  speed  of  sound  as  derived  before  is  given  by 
fx       fx 


C  =  * 


n  ±  4>/360 


In  order  to  define  the  accuracy  with  which  c  is  given  we  use  the 

following  approach: 

log  c  =  log  f  +  log  x  -  log  (n  ±  4>/360) 
* 

d£    =  dJl+dx..  dpi  ±  4./360) 
C        f    x      n  ±  ^/360 

then  the  maximum  error  is  given  by 

M.  =     M.  +  A*        A(n  ±  <j)/360 

c  f        x  n  ±  ^/360 

This  error  will  be  calculated  for  the  highest  (80  kHz)  and  the  lowest 

(30  kHz)  frequencies.  The  coherent  decade  frequency  synthesizer 

synthesizes  the  frequency  of  the  output  signal  by  proper  combination 

of  a  number  of  frequencies,  each  derived  on  a  proportional  basis  from 

a  single  internal  master  frequency  crystal.  It  was  essential  in  this 

application  because  of  the  requirement  for  a  stable,  sine-wave  signal 

source  capable  of  ultraprecise  frequency  adjustments  through  a  wide 

frequency  range.  Since  the  output  is  always  coherent  in  frequency  with 

that  of  the  quartz-crystal  oscillator  source  it  has  the  same  percentage 

accuracy  as  the  source  which  is  essentially  the  smallest  digit  in  the 

specific  frequency  or  1.6  x  10   at  60  kHz.  This  then  implies  that 

~2 
Af     =    10  ■=  =  3.3  x  10"7  at  30  kHz   and 
T       30x1 06 


_2 

^t    =  ]-^-T-    =  1.25  xl0~7    at  80  kHz 

T        80x1  (T 

Next  the  ax/x  term  will  be  examined.  The  absolute  hydrophone  distance 

3 
could  only  be  measured  to  one  part  in  1.5  x  10  (0.5  mm)  which  is 

equivalent  to 

&  =     Jx]0l_      __     64xl0-4 

x  .78 

Assuming  the  Af/f  and  the  A^/^  terms  to  be  zero  this  error  in  measuring 
the  separation  between  the  hydrophones  results  in  an  absolute  error  for 
the  speed  of  sound  of  1  m/sec  (a  level  shift  of  1  m/sec  in  the  non- 
bubbly  value  of  c),  but  does  not  affect  the  dispersion  curve,  that  is 
the  pattern  of  the  speed  of  sound  versus  frequency.  Therefore  this 
error  is  no  serious  problem  and  will  be  ignored  in  the  analysis  since 
we  are  not  interested  in  the  absolute  accuracy  but  rather  the  dispersion 
in  the  speed  of  sound. 

To  reduce  possible  errors  due  to  changes  in  the  geometrical  config- 
uration, especially  in  the  separation  distance  between  the  hydrophones, 
it  was  necessary  to  use  as  a  mounting  a  wire  put  under  100  lb  of  tension 
The  remaining  possible  vibration  due  to  forces  at  sea  was  estimated  to 
be  0.1  mm  and  results  in 

^         -    ifil  =  1.28  xKT4 
x       ,78 

The  accuracy  of  the  phase  was  the  limiting  factor.  Since  basically 
a  phase  interference  approach  was  used  to  measure  the  speed  of  sound 
any  relative  phase  changes  could  not  be  allowed  to  occur  in  the  two  arms 
of  the  system.  It  is  a  question  whether  one  can  assume  the  two  parallel 


parts  within  the  system  (hydrophone-preamplifier-filter-input  channel 
phase  meter)  to  be  identical  with  respect  to  their  influence  on  the 
signal.  Actually,  complete  identity  can  never  be  obtained  over  a  wide 
frequency  range.  Therefore  special  attention  was  given  to  the  phase 
characteristic  in  selecting  the  hydrophones.  The  LC-5  resp. 
LC-10  were  selected  because  their  resonance  is  far  above  the  frequency 
range  of  interest.  Comparing  the  dispersion  curve  for  4.3  m  with  that 
at  14.3  m  also  points  out  that  there  is  no  phase  variation  in  the 
pattern  that  shows  up  consistently  at  all  depths  and  may  therefore  be 
due  to  a  different  phase  behaviour  of  the  two  hydrophones.  In  addition 
the  zero  control  of  the  phase  meter  has  been  used  to  compensate  at  a 
single  frequency  for  unequal  phase  in  the  two  receiving  parts  when  the 
two  input  signals  were  known  to  be  exactly  in  phase  by  using  the 
electrical  signal  from  the  output  of  the  frequency  synthesizer. 

The  accuracy  of  the  analog  output  of  the  phase  meter  is  ±0.1°±0.3% 
of  the  phase  angle  measured  for  all  frequencies  or  ±0.16°  for  the  20° 
phase  angle  which  was  used  normally  in  the  experiment.  This  causes  a 
fractional  error 

A(n±  */360)  =  JS =  .272  x  1(T4  at  30  kHz 

n  ±  4>/360     (360  x  16)  +  20 

.108  x  10~4  at  80  kHz 


(360  x  41)  +  20 


There  are  of  course  accuracy  limitations  in  measuring  the  phase  of 
distorted  wave  forms.  Because  the  phase  meter  makes  its  phase  measure- 
ments by  measuring  the  fraction  of  a  cycle  between  the  negative-going 
zero  axis  crossing  of  the  two  test  signals  harmonic  or  scattered  signals 


can  displace  the  zero  crossing  in  either  a  plus  or  a  minus  manner 
depending  on  the  phase  relationship  to  the  test  signal. 

The  scattered  signal  was  measured  to  be  down  by  37.4  dB  re  1  V/y bar 
as  compared  with  the  main  signal  which  for  the  worst  case  could  result 
in  an  error  of  .77°.  A  scattered  signal  a'  adding  to  the  main  signal  § 
would  result  in  the  largest  error  in  e  when  the  condition  depicted 
below  occurs 


such  that         20  log  x-  =  37.4  dB 

a 


for  small   e  sine^6     =     l- 


then 


b 

b-    =     lO37'4/20     -     74 
a 


and  e     =     jj    =     .0135  in  radians 

2  ,77° 


A  possible  source  of  error  is  also  the  fact  that  at  some  frequencies 

the  response  of  one  hydrophone  is  significantly  stronger  than  the  other 

which  means  the  possibility  of  harmonic  content  in  the  signal  at 

multiples  or  submultiples  of  these  frequencies.  The  limit  error  for 

a  harmonic  of  magnitude,  e,  is  0  =  —  57°  where  a  is  the  magnitude  of 

a 

the  fundamental  (manual  for  digital  input  phase  meter  model  355). 
Since  the  second  harmonic  was  measured  to  be  -46.34  dB  re  1  V/ybar 
as  compared  with  the  fundamental  the  harmonic  content  could  at  most 
result  in  an  error  of  0.28°. 

Assuming  the  worst  condition  for  the  phase  reading  or  an  error  of 
.16°  at  the  phase  meter,  a  maximum  additional  error  of  .77°  due  to  the 


scattered  signal,  and  also  a  maximum  error  of   .28°  due  to  harmonic 
content  in  the  signal   the  total   phase  error  could  be 

A(n  ±  ^360)     -  J6+   -77+   -28     ■     2.08  x  10-4  at  30  kHz 

n  ±  <f>/360  (360  x  16)  +  20 


1  207  -4 

.82  x  10  *  at  80  kHz 


(360  x  41 )  +  20 
Summarizing  it  can  be  stated 


—  =  3.3  x  10~7  +  1,28  x  10~4  +  2.08  x  10"4  =  3.36  x  10"4  at  30  kHz 
c 

=  1.25-xlO'7  +  1.28  x  10~4  +  .82  x  10"4  =  2.1  x  10~4  at  80  kHz 


Therefore  the  maximum  error  in  c  is  mainly  determined  by  the  percentage 
error  in  x  and  in  <j>  and  is  given  by 

Ac  =  (1500)  (3.36  x  10~4)  =  .5  m/sec  at  30  kHz  and 

=  (1500)  (  2.1  x  10"4)  =  .31  m/sec  at  80  kHz 

Comparing  the  value  for  Ac  at  the  lower  and  at  the  upper  limit  points 
out  that  in  general  one  is  able  to  take  more  accurate  measurements  at 
the  higher  frequencies. 


XI.  SUMMARY  AND  CONCLUSIONS 

A  fixed  source-hydrophone  system  separated  by  1.8  m  was  used  to 
study  a)  the  in-situ  speed  of  sound  as  a  function  of  frequency  and 
depth  and  b)  the  temporal  characteristics  of  the  phase  fluctuations 
of  a  continuous  wave  acoustical  signal  over  the  frequency  range  25  -  80 
kHz.  The  experiment  was  carried  out  during  Sea  State  1  conditions  at 
the  NUC  Oceanographic  Research  Tower  in  water  of  depth  56  ft  one  mile 
off  shore  at  Mission  Bay/San  Diego  in  October  1971.  The  near  surface 
region  is  of  particular  interest  because  bubbles  and  turbulence  can 
have  significant  acoustical  effects  in  addition  to  those  caused  by 
temperature  and  salinity  inhomogeneities. 

Four  different  depths,  4.3,  7.3,  9.3,  and  14.3  m,  were  examined. 
The  near  surface  values  of  c  ranged  from  +  6  m/sec  to  -3  m/sec  relative 
to  the  value  for  bubble-free  sea  water,  with  a  maximum  error  of  0.5  m/sec 
The  speed  of  sound  was  found  to  peak  at  approximately  35  and  71  kHz  at 
4.3  m  depth.  The  lower  peak  frequency  increased  as  the  depth  is  in- 
creased but  not  at  the  rate  which  would  be  expected  for  bubbles  that 
are  wind-initiated  at  the  surface  and  carried  to  greater  depths  without 
losing  their  identity.  A  possible  explanation  is  that  the  bubbles 
originated  within  the  volume  or  at  the  bottom  and  lost  gas  by  diffusion 
as  they  rose  to  the  surface.  Another  possibility  is  that  presence  of 
more  high  frequency  bubbles  at  the  surface  affects  the  apparent 
resonance  frequency  of  the  larger  bubbles  as  determined  by  sound  dis- 
persion. The  gradual  and  smooth  drop-off  in  the  resonant  behaviour  with 
increasing  depth  demonstrates  that  bubbles  of  this  average  size  are 
somewhat  continuously  distributed  in  the  medium  down  to  14.3  m.  The 


second  distinct  peak  centered  at  71  kHz  appears  to  be  only  a  character- 
istic of  the  shallowest  run  and  does  not  repeat  at  the  other  depths. 
It  is  postulated  that  these  small  bubbles  are  introduced  at  the  surface 
either  by  randomly  breaking  waves  or  by  continental  aerosols  that  drop 
into  the  water  and  carry  tiny  air  bubbles  with  them  (Medwin,  1970). 
The  surface  initiated  bubbles  disappear  quite  fast  due  to  gas  diffusion 
out  of  the  bubble  when  transported  downwards  in  the  medium.  We  conclude 
from  the  foregoing  that  we  encountered  bubbles  of  radii  centered  around 
124  and  54  microns  and  resonating  at  30  and  69  kHz  respectively.  The 
broad  dispersive  characteristic  around  30  kHz  implies  many  bubbles 
larger  than  60  microns  whose  principal  source  was  at  the  surface  and 
which  were  entrained  by  convective  currents  carrying  them  downwards. 

Statistical  analysis  shows  that  for  frequencies  other  than  at  the 
peak  bubble  populations  the  phase  fluctuations  can  be  approximated  by 
a  Gaussian  probability  density  function.  However  around  the  peak 
variations  in  sound  speed  the  PDF  loses  the  smooth  Gaussian  appearance 
and  shows  a  wide,  somewhat  irregular  spread.  The  temporal  variation 
of  the  standard  deviation  of  the  sound  phase  was  small  (i  .32°).  How- 
ever, the  standard  deviation  of  the  phase  fluctuation  changed 
considerably  as  a  function  of  frequency  increasing  to  1.92°  at  40.6 
kHz  and  to  2.56°  at  70  kHz,  the  two  sound  dispersion  centers  of  this 
experiment. 

Typically,  the  phase  was  decorrelated  O/e)  after  1.46  sec.  This 
compares  with  preliminary  oceanographic  information  which  shows  corre- 
lation times  of  3,96  sec  for  the  salinity  and  2.28  sec  for  the  speed 
of  sound  as  measured  by  the  velocimeter.  These  values  show  internal 
consistency. 


The  power  spectral  density  of  the  varying  phase  shows  its  strongest 
values  at  frequencies  less  than  0.5  Hz;  this  is  presumably  due  to 
bubbles  entrained  at  orbital  frequencies  associated  with  the  surface 
wave  system  which  has  maximum  energy  in  the  same  frequency  range. 

Due  to  time  considerations  it  was  not  possible  to  further  examine 
the  phase  fluctuations  in  terms  of  the  functional  dependence  on 
temperature,  salinity,  turbulence,  and  wave  height.  The  analysis  of 
the  oceanographic  parameters  will  be  carried  on  in  the  theses  by 
Lt.  Bordy,  LCDR  Duchock,  and  Lt.  Seymour  (March  1972)  and  the  emphasis 
should  then  be  placed  on  the  interrelation  of  the  parameters  as  given 
by  the  cross-correlations  or  cross  spectra  for  instance. 


XII,  RECOMMENDATIONS 

It  is  recommended  that  research  in  this  area  be  continued,  especially 
by  making  an  in-situ  determination  of  the  bubble  concentrations  at  the 
same  time  as  sound  speed  measurements  are  taken  in  order  to  determine 
the  variation  of  bubble  population  with  location,  weather  conditions, 
day  or  night,  and  local  effects.  It  is  also  suggested  to  examine  a 
smaller  frequency  range  in  more  detail  rather  than  to  skim  through  a 
wide  frequency  range  in  larger  steps  and  thereby  losing  valuable 
information.  Frequencies  below  30  kHz  seem  to  be  most  inviting  because 
of  relevance  to  the  Navy  and  promising  because  the  largest  changes  in 
the  speed  of  sound  occur  here,  just  below  the  frequency  of  the  most 
prominent  bubble  population. 


APPENDIX  A.  PRESENTATION  OF  DATA 


Run  6, 

4.3  m 

freq 
[kHz] 

real   t 
[sec] 

mean 
[Dig.] 

standard 

dev. 

[Dig.] 

mean 
[radians] 

n 

c=^  x=,783 

30.72 

180 

-.84 

.06 

-.07467 

15.0987 

1511  .99 

32.70 

101 

-.83 

.08 

-.07378 

16.9096 

1514.18 

34.71 

90 

-.69 

.06 

-.06133 

17.9220 

1516.46 

36.82 

80 

+  .445 

.11 

+.03956 

19.0562 

1512.90 

38.70 

30 

+  .65 

.08 

+.05778 

20.0745 

1509.48 

40.60 

38 

+  .83 

.12 

+.07378 

21.0905 

1507.30 

42.42 

44 

+  .08 

.08 

+.00711 

22.0238 

1508.13 

46.40 

193 

+  .44 

.06 

+.03911 

24.0558 

1510.29 

48.40 

107 

+  .665 

.06 

+.05911 

25.0758 

1511  .31 

50.30 

102 

+  .48 

.06 

+.04267 

26.0593 

1511.36 

52.22 

96 

+  .485 

.03 

+.04311 

27.0598 

1511.03 

54.14 

102 

+  .49 

.06 

+.04356 

28.0602 

1510.74 

56.09 

99 

+  .57 

,06 

+.05067 

29.0673 

1510.92 

58.00 

98 

+  ,49 

.06 

+.04356 

30.0602 

1510.77 

59.92 

90 

+  .54 

.06 

+.04800 

31.0647 

1510.31 

61.82 

89 

+  .53 

.06 

+.04711 

32.0638 

1509.65 

63.72 

89 

+  .50 

.07 

+  ..  04444 

33.0611 

1509.11 

65.64 

103 

+  .51 

,11 

+.04533 

34.0620 

1508.90 

67.59 

121 

+  .58 

,11 

+.05156 

35.0682 

1509.14 

69.60 

104 

+  .66 

.14 

+,05867 

36.0753 

1510.64 

71.56 

86 

+  .65 

..10 

+.05778 

37.0745 

1511.32 

73.49 

90 

+  .70 

,11 

+,06222 

38.0789 

1511.14 

75,38 

101 

+  .59 

.09 

+.05244 

39.0691 

1510.72 

77,28 

110 

+  .56 

.07 

+.04978 

40.0665 

1510.25 

79.20 

78 

+  .57 

.07 

+,05067 

41,0673 

1510.05 

Run  7, 

4.3  m 

freq 
[kHz] 

real   t 
[sec] 

mean 
[Dig.] 

standard 

dev. 

[Dig.] 

mean 
[radians] 

n 

fx 
C=~  x=.783 
n 

30.69 

73.6 

-.94 

.06 

-.0836 

15.8997 

1511.37 

32.69 

112.8 

-.85 

.04 

-.0756 

16.9077 

1513.88 

34.69 

113.6 

-.87 

.05 

-.0773 

17.9060 

1516.94 

36.60 

88.0 

-.97 

.04 

-.0862 

18.8971 

1516.52 

38.76 

98.4 

+  .65 

.09 

+.0578 

20.0745 

1511.82 

40.61 

93.6 

+  .43 

.07 

+.0382 

21.0549 

1510.22 

40.61 

93.6 

+  .44 

,07 

+.0391 

21.0558 

1510.16 

42.53 

104.8 

+  .53 

.05 

+.0471 

22.0638 

1509.30 

44.45 

93.6 

+  .47 

.04 

+  .0418 

23.0585 

1509.39 

48.34 

88.8 

+  .55 

.05 

+.0489 

25.0656 

1510.05 

52.19 

92.8 

+  .56 

.06 

+.0498 

27.2667 

1509.79 

56.04 

104.8 

+  .48 

.06 

+.0427 

29.0594 

1509.99 

59.86 

78.4 

+  .45 

.09 

+.0400 

31.0567 

1509.19 

63.66 

99.2 

+  .51 

.09 

+.0453 

33.0620 

1507.65 

67.50 

99.2 

+  .47 

..07 

+.0418 

35.0585 

1507.55 

71,44 

73.6 

+  .52 

.09 

+.0462 

37.0629 

1509.26 

75.24 

92.0 

+  .50 

.08 

+.0444 

39.0611 

1508.22 

79.10 

105.6 

+  .51 

.07 

+.0453 

41.0620 

1508.34 

Run  8, 

9.3  m 

freq 
[kHz] 

real   t 
[sec] 

mean 
[Dig.] 

standard 

dev 

[Dig.] 

mean 
[radians] 

n 

c=%  x=.783 
n 

25.17 

81.6 

+  .51 

.03 

+.0453 

13.062 

1508.81 

28.98 

79.2 

+  .60 

.03 

+.0533 

15.070 

1505.73 

32.91 

84.8 

+  .56 

.05 

+.0498 

17.0665 

1509.89 

36.79 

96,0 

+  .46 

.03 

+.0409 

19.0576 

1511.55 

40.62 

91.2 

+  .54 

,03 

+.0480 

21.0647 

1509.89 

44.49 

93.6 

+  .57 

.06 

+.0507 

23.0674 

1510.17 

48.34 

102.4 

+  .59 

.07 

+.0524 

25.0691 

1509.84 

52,18 

95.2 

+  .59 

.04 

+.0524 

27.0691 

1509.36 

56.07 

103.2 

+  .51 

,10 

+.0453 

29.0620 

1510.66 

59.94 

88.8 

+  .58 

.08 

+.0516 

31.0683 

1510.64 

63.71 

101.6 

+  .48 

.08 

+.0427 

33.0594 

1508.95 

67.56 

108.0 

+  .49 

.07 

+.0436 

35.0603 

1507.92 

71.39 

54.4 

+  .46 

.05 

+.0409 

37.0576 

1508.42 

Run  9, 

14.3  m 

• 

freq 
[kHz] 

real   t 
[sec] 

mean 
[Dig.] 

standard 

dev. 

[Dig,] 

mean 
[radians] 

n 

c=^-x=,783 

28.95 

42.4 

+  .70 

.05 

+.0622 

15.0789 

1503.28 

36.76 

28.0 

+  .62 

.08 

+.0551 

19.0718 

1509.20 

44.48 

32.0 

+  .56 

.04 

+.0498 

23.0665 

1509.89 

55.93 

36.8 

+  .58 

.06 

+  .0516 

29.0683 

1506.56 

59,78 

36.8 

+  .54 

.05 

+.0480 

31.0647 

1506.78 

63.57 

25.6 

+  .53 

,06 

+.0471 

33.0638 

1505.43 

71,28 

78.4 

+  .76 

.09 

+.0676 

37.0843 

1505.01 

Run  10, 

7.3  m 

freq 
[kHz] 

real   t 
[sec] 

mean 
[Dig.] 

standard 

dev. 

[Dig.] 

mean 
[radians] 

n 

c=^-x=.783 
n 

27,05 

81.6 

.59 

.04 

.0524 

14.0691 

1505.44 

30.96 

71.2 

.70 

.05 

.0622 

16.0789 

1507.67 

34.89 

96,8 

.51 

.04 

.0453 

18.0620 

1512.51 

38,72 

89.6 

.55 

.04 

.0489 

20.0656 

1510.93 

42,55 

87.2 

.61 

.04 

,0542 

22.0709 

1509.53 

46.40 

89.6 

.56 

,05 

.0498 

24.0665 

1509.62 

50.23 

81.6 

.57 

,09 

.0507 

26.0674 

1508.78 

54,09 

86,4 

.52 

,05 

.0462 

28.0629 

1509.20 

57.70 

84,8 

-1,08 

.07 

-.0960 

29.8873 

1511.65 

61.59 

80.0 

-.79 

.06 

-.0702 

31.9131 

1511.13 

65,59 

78,4 

,48 

,13 

.0427 

34.0594 

1507.86 

C      PROGRAM  1 

C      PROGRAM  CONVERK  INPUT, OUTPUT, PUNCH) 

C      AVGX  IS  THE  AVERAGE  DISTANCE  DISTANCE  THE  RECORDING 

C      DRUM  ADVANCES  PER  ONE  HOUR  REAL  TIME,  I.E.  THE  DATA 

C      SAMPLING  RATE, 

DIMENSION  U(6D00), V< 6000 ) , N( 80 ) , IBUFF( 5000 ) ,NK( 80 ) 

IK=0 

READ  97, L 

97  FORMAT ( 12) 
READ  98, AVGX 

98  FORMAT(FIO.O) 
READ  96,IFILE 

96  F0RMAT(I2) 

DO  116  JJ=1,L 
DELT=36.0/AVGX 
PRINT  99,DELT 

99  F0RMAT(1H0,2X,25H  SAMPLING  INTERVAL  EQUALS,  1X,E15.7, 
11X,7HSEC0NDS) 

PRINT  700, JJ 
700  F0RMAT(1H0,9HJJ  EQUALS,  2X , I  2 ) 
DO  100  1=1,5500 
U(I)=0.0 

100  V<I)=0.0 

101  DO  202  1=1,5000 
202  IBUFF(I)=0 

COUNTX=0.0 

COUNTY=0.0 

NUM=5000 

M=l 

K8=0 

CALL  LI0F(5LRBCD1, IBUFF, NUM, NPAR, NEOF ) 

IF(NEOF)  602,200 

200  K=-7 
KF  =  0 

201  K=K+8 
KA=K+8 
KC=0 

DO  102  KB=K,KA 
KC=KC+ 1 
NK(KC)=IBUFF(KB) 

102  IF  (NK(KC) .EQ.O)  KF=1 

DECODE (80, 103, NK)  (N(I),I=1,80) 

103  F0RMAT(80R1) 
DO  104  19=1,30 

IF(N( 19) .EQ.50B)  GO  TO  106 

IFCN(I9).EQ.55B)  GOTO  107 

IF(NCI9).EQ.34B)  GO  TO  108 
C      SYMBOL  /  (50)  REPRESENTS  AN  INCREMENT  TRAVEL  IN  THE 
C      MINUS  X  OR  Y  DIRECTION  BY  THE  DIGITIZER. 
C      SYMBOL  0,(55B),  REPRESENTS  A  ZERO  INCREMENT  TRAVEL 
C      IN  THE  X  OR  Y  DIRECTION  BY  THE  DIGITIZE11- 
C      SYMBOL  1,(345),  REPRESENTS  AND  INCREMENT  TRAVEL  IN 
C      THE  POSITIVE  X  OR  Y  DIRECTION  BY  THE  DIGITIZER 
C      SYMBOL  *v(47B)f  IS  A  FLAG  INSERTED  IN  IN  THE  RECORD 
C      BY  THE  PERSON  DIGITIZING  BY  USE  OF  THE  IRG 

GO  TO  104 

106  RX=-0.01 
K8=K8+1 
GO  TO  109 

107  RX=0.0 
K8=K8+1 
GO  TO  109 

108  RX=0.01 
K8=K8+1 

109  K3=K8/2 

1/    T  =  O   J;  1/   -2 

IF(K3.EQ.K8)  GO  TO  111 
COUNTX=COUNTX+RX 
IF(COUNTX.'MF.O.O)  GO  TO  110 
GO  TO  104 

110  U(M)=COUNTX 
GO  TO  104 


111  COUNTY=COUNTY+RX 
IF(COUNTX.NE.O.O)  GO  TO  112 
GO  TO  104 

112  V(M)=COUNTY 
COUNTX=0.0 
M=M+1 
IF(M.GT.6000)  GO  TO  600 

104  CONTINUE 
IF(KF.EQ.l)  GO  TO  113 
GO  TO  201 

113  MAX=M-7 

105  PRINT  205, M, 19 
205  FORMAT( 1H0,2I10) 

C      TOTAL  TIME  OF  THE  RECORD  EQUALS  NUMBER  OF  DATA 
C      POINTS  TIMES  THE  SAMPLING  INTERVAL,  DELT. 

TIME=(M*DELT)/ 3600.0 

PRINT  115, TIME 

115  FORMAT (1H, 10X,28H  TOTAL  TIME  OF  RECORD  EQUALS, 2X, 
1E15.7,2X,7H  HOURS. ) 

PRINT  215, (V(I ) ,1=1, M) 

PRINT  215, (U(I ),I=1,M» 
215  FORMATQH  ,10X,14F7.2) 

PUNCH  213,(V(I),I=1,M) 
213  F0RMATU4F5.2) 

PUNCH  2100 
2100  FORMAT  ( 20H******************** ) 

116  CONTINUE 
STOP 

600  PRINT  601 

601  F0RMAT(1H0,20X,37H*****  U  AND  V  SPACE  INADEQUATE  *** , ) 

602  IK=IK+1 
IF(IK.GT.IFILE)603,200 

603  STOP 
END 


C      PROGRAM  2 

C      PROGRAM  FOR  THE  SPECTRAL  ANALYSIS  OF  APERIODIC  RECORDS 

C      DATA  MUST  BE  SAMPLED  AT  EQUAL  INTERVALS  OF  DT 

C      CALCULATES  EITHER  THE  AUTO  SPECTRA  OR  THE  CROSSSPECTRA 

C      IF  NFLAG1  =  0  ,  AUTOSPECTRA IF  NFLAG1=1     COSPECTRA 

C      IF  NFLAGP  =  1  CALL  PRESS,  TOTE 

C       IF  NFLAGC  =  1  CALL  URMS 

C  NTS=  NUMBER  OF  TIME  SAMPLES 

C  MLAG  =  NUMBER  OF  TIME  LAGS 

C  DT  =  TIME  INCREMENT 

C  DHZ  =  FREQUECY  SPACING 

C  FBHZ  =  LOWEST  FREQUENCY  (HZ) 

C  FEHZ  =  HIGhEST  FREQUENCY  OF  INTEREST  (HZ) 

C  TM  =DT*MLAG 

C  DHZ  =  1/(2*A*TM) 

C  FEHZ  =  B/(2*TM)=  B/ ( 2*MLAG*DT ) 

C  FEHZ  =  A*B*DHZ 

C  A*B  =  TOTAL  NO.  FREQUENCY  SPACINGS  OF  ENERGY  SPECTRUM 

C  A  AND  B  ARE  INTEGER  CONSTANTS 

C  TRANS —  INSTRUMENT  AND  REMARKS 

REAL*8  ITITLE(12),ITITEL(12),ITITIL(12) 

DIMENSION  cl(6000) ,F2(6000) 

DIMENSION  PHI (600) ,TAU(600) ,SPHI (600) , APHI ( 600) , 
1PHN(600) 

DIMENSION  FREQ(700) ,CYCL ( 700 ), SPEC ( 700) ,SSP(700) 

DIMENSION  CSPEC(700),QSPEC(700) ,SS2(700) 

DIMENSION  RESPF(700),PER(700),SPE2(700) , PHASE (700), 
1C0HER(700) 

DIMENSION  DATE (2), HOUR (2), TRANS (10) 

REAL*4  LABEL/4H  CO  /, LASLI /4HQUAD/, MONT , LABI L/4H  PHI/, 
1LABLE/4HCR0S/ 

READ(5,200,END=8  88)  ( I TITLE! I ) , 1=1 ,12 ) 

READ(5,2  00,END=8  88)     (ITITEL(I)  ,1  =  1,12) 
200    F0RMAT(6A8) 

READ(5,802)  DATE , HOUR, TRANS 

802  F0RMAT(2A4,2A4,10A4) 
READ(5,96)  NFL AG1 , NFL  AGP, NFL AGC 
READ(5,99)  NTS , MLAG, DT , DHZ , FBHZ , FEHZ 
RE AD (5, 8  01)  CALX1,CALX2,H,X2,AZZ,S 
WRITE(6,811) 

811  FORMAT(lHl) 

WRITE(6,803)  DATE, HOUR , TRANS 

WRITE (6, 98)  NTS, ML AG, DT, DHZ, FBHZ, FEHZ, CALX 1,CALX2,H,X2 
96  FORMAT  (316) 

98  F0RMAT(7H  NTS  =  ,16//, 8H  MLAG  =  ,I6//,6H  DT  =  ,F9.5//, 
16H  DHZ=  ,F10 

2.8//,  ■  FBHZ  =  SFIO.S//,1  FEHZ  =  ,,F12.8//,*  CALX1  =  ■ 

3,F10.8//, 

4»  CALX2  =  ',F10.8//,«  H  =  «,F10.3//,»  X2  =  • .F10.3//I 

99  FORMAT  (2I6,5F10.8) 
801  FORMAT(6F10.5) 

803  F0RMAT(40H  SOUND  PHASE  STUDY  OF  SAN  DIEGO  BAY      , 
15X,7H  D*TE   ,2A4,10H    HOUR    , 2 A4// , 5X, 1 0A4/// ) 

C  BAND  WIDTH  FREQUENCIES  OF  TOTAL  ENERGY  FLUX-  CMIN,CMAX 

CMIN  =  0.0 

CMAX  =  0.2 
C 
C  COMPUTING  POWER  SPECTRUM 

PI=3. 14159265 

FB  =  FBHZ-2.0*PI 

FE  =  FEHZ-2.0-PI 

OF  =  DHZ  *2.0*PI 

FMIN  =  2.0*PI*CMIN 

FMAX  =  2.0*PI*CMAX 
900  FORMAT( 14F5.2) 

READ(5,900)(F1( I), 1=1, NTS) 

XMAX    =    Fl(l) 

DO    23    1=2, NTS 

IFCXMAX.GE.FK I ) )    GO    TO    23 

XMAX  =  FKI) 
23  CONTINUE 

DO  24  1=1, NTS 


F1(I)  =  (F1(  D-XMAX) 
24    CONTINUE 

WRITE(6,901)(F1( I), 1=1, NTS) 
901    F0RMAT(14F7.2) 

CALL    TREND(F1 ,NTS, DT,C ALX1 ) 
21    DO    10    M=1,MLAG 

SUM=0.0 

NMAX=NTS-M+1 

DO  8  I=1,NMAX 

NN=M+I-1 
8  SUM=SUM+F1( I)*F1(NN) 

XNMAX=NMAX 

XX=M-1 

TAU(M)=XX*DT 

PHI(M)=SUM/XNMAX 

PHN(M)     =    PHI(M)/PHI(1) 
10    CONTINUE 

CALL    PARZ(MLAG»PHI ) 
798    FORMATC     PHI(O)     =    'tFlO.4//) 

NFREQ=(FE-FB)/DF+0.1 

DO    14    N=1,NFREQ 

XN=N 

FREQ(N)=(XN-1.0)*DF+FB 

14  CYCL(N)    =    FREQ(N)/(2.0*PI) 
PER(l)     =    0.0 

DO    15    N    =    2,NFREQ 

15  PER(N)     =    1.0/CYCL(N) 
MLAGM1=MLAG-1 
XMLAG=MLAG 

46    DO    50    N=1,NFREQ 

SUM=0.5*(PHI(1)+PHI(MLAG)*C0S(FREQ<N)*TAU(MLAG)  )  ) 

C1=C0S(FREQ(N)*DT) 

S1=SIN(FREQ(N)*DT) 

CC=1.0 

SC=0.0 

DO  49  M=2,MLAGM1 

CT=CC*C1-SC*S1 

ST=SC*C1+CC*S1 

CC=CT 

SC=ST 

49  SUM=SUM+PHI (M) *CC 

50  SPECCN»=SUM*2.0/XMLAG 
IF(NFLAGP)     741,741,742 

742  CONTINUE 

CALL  PRESS (FREOf SP EC , NFREQ,H , X2,DF, FB ,FMAX) 
741  CONTINUE 

MD=3.1415926  5/(DF*DT*XMLAG)+0.1 

DC  =  DF/(2.0*PI ) 

SUM  =  0.0 

DO  650  N=1,NFREQ 

650  SUM  =  SUM  +  SPEC(N) 
WRITE(6,651)  SUM 

651  FORMAT  { 25H      VARIANCE  OF  SPECTRUM//, 3X, 
116H      VARIANCE  =  ,F10.5,10H    M2      //) 

WRITE(6,106) ( TAU(M) , PHN(M) ,M  =  1 ,MLAG) 
INUM=FEHZ/DHZ 
DO  9988  1=1,  I  MUM 
SPEC ( I )=10*AL0G10(SPEC(I ) ) 
9988  CONTINUE 
C      CALL  DRAW  TO  PLOT  SPECTRUM  OF  SAN  DIEGO 
C      AS  A  FUNCTION  OF  FREQUENCY 

CALL  DRAW(200,TAU,   PHN, 0, 0, LABEL , I T ITEL , 8 . 0 ,0.0, 2, 0, 2 
1,0, 8, 8,1, LAST) 

CALL  DRAW( I NUM,CYCL , SPEC ,0 ,0 , LABEL, IT ITL E , 0 .0 ,0. 0 ,0 , 0, 
10,0,8,8,1, LAST) 

105  FORMAT  (10X,2^H  ENERGY  SPECTRAL  DENS ITY/// , 5X, 
110H  FREQUENCY, 5X,13H 

2  RAW  SPECTRUM, 3X,18H  SMOOTHED  SPECTRUM//) 

106  FORMAT  ( 5X , F10 . 5 ,6X ,F10. 5 ) 
WRITE(6,105) 

WRITE(6, 106)(CYCL( M),SPEC(M)        ,M=1,NFREQ) 
888  RETURN 


END 

SUBROUTINE  PRESS ( FREQ, SPEC, NFREQ, H,X2,DF , FB, FMAX) 
C  SUBROUTINE  TO  CONVERT  PRESSURE  SPECTRUM  TO  ENERGY  SPECTRUM 

DIMENSION  FREQ( NFREQ), SPEC (NFREQ) 

NMAX  =  (FMAX-FBJ/DF+1.0 

DO  11  1=1, NMAX 
C  CALCULATE  LINEAR  WAVE  LENGTH  BY  NEWTONS  METHOD 

IF(FREQ( I)-0. OOOOl)  8,8,7 

7  XKHO  =  FREQ(I)*FREQ(I)*H/9.80 
IFtXKHO-6.3)  5,1,1 

1  XKH  =  XKHO 

GO  TO  9 
5  XKH  =  SQRT(XKHO) 

3  SH  =  SINH(XKH) 
CH  =  COSH(XKH) 

EPS  =  XKHO-XKH*SH/CH 
SLOPE  =  -XKH/CH*-2-SH/CH 
DXKH  =  -EPS/SLOPE 
IF(ABS(DXKH/XKH)-0.0001)  9,9,4 

4  XKH  =  XKH+DXKH 
GO  TO  3 

8  RESPF  =  1.00 
GO  TO  11 

9  XK  =  XKH/H 

RESPF  =  C0SH(XK*H)/(C0SH(XK*(H-X2))) 
IF(RESPF-IO.O)  11,11,12 

11  SPEC(I)  =  RESPF*RESPF*SPEC(I) 

12  RETURN 
END 

SUBROUTINE  TREND ( FX,NTS,DT,CALXX) 

DIMENSION  FX(NTS) 
C  CALIBRATIONG  RECORD 

DO  104  1=1, NTS 
104  FXU)  =  FX(  I)*CALXX 
C  COMPUTING  THE  LINEAR  TREND 

FNTS  =  NTS 

SUMF  =  0.0 

DO  101  1=1, NTS 

101  SUMF  =  SUMF  +  FXU  ) 
SUMF1  =  0.0 

DO  102  1=1, NTS 
XI  =  I 

102  SUMF1  =  SUMF1  +  XI*FX(I) 
XNM1  =  NTS-1 

XNP1  =  NTS+1 

XM  =  (1.0/DT)*(12.0*SUMF1/(FNTS*XNM1*XNP1)-6.0*SUMF/ 

1(XNM1*FNTS) ) 
B   =    SUMF/FNTS-XM*XNPl*DT/2.0 
FMEAN    =     SU'«1F/FNTS 
WRITE(6,9)     FMEAN, XM,B 
9    FORMAT     (3X,8H    MEAN    =    ,F10.5,3X,9H    SLOPE    =    ,F10.5,3X, 
113H    INTERCEPT=     ,F10.5//) 
DO    103    1=1, NTS 
XI    =    I 

103  FXU)    =    FXU)    -    (B  +  XM*XI*DT) 
RETURN 

END 

SUBROUTINE    SMO ( MD, XI , X2, NFREQ ) 

DIMENSION    XKMD)  ,X2(MD) 

DO    1    N=1,MD 

NA=N+MD 

NN=NFREQ-N+1 

NB=NN-MD 

X2(N)  =  0.25*(X1(1)+X1(NA) )+0.5*Xl(N) 
1  X2(NN)=0.5*tXl(NN)+XKNB)J 
3  MB=MD+1 


ME=NN-1 
5  DO  2  N=MB,ME 

NA=N+MD 

NB=N-MD 
2  X2(N)=0.25*(X1(NA)+X1(NB))+0.5*X1(N) 

RETURN 

END 

SUBROUTINE  PARZ ( ML AG, PHI ) 

PARZ  SUBROUTINE  PARZEN  FILTERS  AUTO-CORRELLATION  FUNCT 

DIMENSION  PHKMLAG) 

XMLAG  =  MLAG 

MLAGH  =  XMLAG/2. O-O.l 

MLAGH1  =  MLAGH  +  1 

DO  31  M=l, MLAGH 

MM  =  M-l 

R  =  MM 

RM  =  R/XMLAG 

UM  =  1.0-6,0*RM*RM*(1.0-RM) 

PHKM  )  =  PHKM  )*UM 

31  CONTINUE 

DO  32  M  =  MLAGH1,MLAG 

MM  =  M-l 

R  =  MM 

RM  =  R/XMLAG 

RM1  =  (  1.0-RM) 

UM  =  2.0*RM1*RM1*RM1 

PHKM  )  =  PHKM  )*UM 

32  CONTINUE 
RETURN 
END 


C     PROGRAM  3  CALCULATES  DISPERSION  IN  C 
//   EXEC  FORTCLGP,  REGION.  GOIOOK 
//FORT.SYSIN  DO  * 

REAL*8  TITLE(12)/' BOX  209  »,U*«         •/ 
REAL  LABEL/'     •/ 

DIMENSION  A(10),B( 10 ) , MR! 10 ) , W ! 500 ) ,CC  {  500 ) 
DATA  PI/3.1415926/ 
C  SET  MAX  DIMENSION  TO  WORK  WITH. 

N0UT=500 
C  READ  PZ,RHZ,CZ. 

READ(5t100)  PZ,RHZ,CZ 
100  F0RMAT!3F10.5) 
PZ=PZ*1QOOOO.O 
C  READ  RADIUS  (MICRONS)  AND  FPEQ. 

READ (5, 200)  IS R, I  OR, IE R,  IS W, IDW, IEW 
200  FORMAT(6I10) 
C  READ  NO.  CF  SEGMENTS  AND  THEIR  VALUEX...MAX  OF  10  SEGMENTS 
READ(5t300)  NS 
300  FORMAT (12) 

READ(5,400)(A(  I ) , B (  I ) , MR( I ) , 1  =  1 , NS ) 
400  FORMAT(2F10.5,I10) 
C  PRINT  INPUT. 

WRITE(6,500) 
500  FORMAT (  ■  1INPUT:  '  ) 

WRITE(6,510)  PZ,RHZ,CZ 
510  FORMAT (•  PO  =  »,E16.8,T28 
1E16.8) 
WRITE(6,520)  ISR,IDR,IEk 
520  FORMAT('  STARTING  RADIUS  = 
1T54, 'ENDING  RADIUS  =  ',16) 
WRITE!6,530)  ISW,IDWfIEW 
530  FORMAK*  STARTING  W  =  ',16 
1» ENDING  W  =  ',16) 
WRITE(o,540) 
540  FORMAT (»    SEGMENT ■ ,7X ,  • A' 
WRITE(6,550)  ( I , A{ I  J ,B I I ) , 
5  50  FORMAT ( 5X, I2,3X , El 6. 8 , 2X,E 16 . 8 ,4Xj I  6) 
C  CHECK  FOR  INPUT  ERROR.  MR  SHOULD  INCREASE. 
NJ=NS-1 

IF  (NJ  .EQ.  0)  GO  TO  7  00 
DO  600  J=1TNJ 

IF  (MRU)  .LT.  MRU  +  I))  GO  TO  600 
WRITE(6,560) 
560  FORMAT!  •  SEGMENT  RADII  ARE  OUT  OF  SEQUENCE.*) 

STOP 
600  CONTINUE 
C  NOW  CALCULATE  SOME  CONSTANTS  TO  SAVE  TIME. 
700  IR=(IER-ISR)/IDR+1 
IW=( IEW-ISWJ/IDW+1 
WC=1.975*SQRT(PZ/RHZ)*1000000.0 
PIC=0.C27/( (PI -2000. 0)**0.322) 
C  CLEAR  OUTPUT  BUFFER 
DO  800  I=1,N0UT 
W(  I  )=0. 
800  CC(I)=0. 
C  INITIALIZE  STARTING  VALUES. 
C  WW=WORKING  OMEGA. 

C  WZI=CALCULATED  OMEGA  ZERO  SUB  I. 
C  RW=WORKING  RADIUS. 
ITEMPW=ISW 
DO  2000  K=1TIW 
WW=ITEMPW 
W(K)=ITEMPW 
IF  (K.NE.l )  GO  TO  915 
WRITE(6,900)  ITEMPW 
900  FORMAT ( • 1W  =  • ,16) 

WRITE<6,910) 
910  FORMAT!'    INDE X • , 6X, « RADI US • , 1 IX , • WZ ■  , 1 3X  ,  ■ U ! R) D < R ) • ) 
915  SUM=0.0 

ITEMPR=ISR 
C  THIS  IS  THE  SUMMATION  LOOP. 
DO  1000  L=1,IR 


•RHO 

=  »,E16.8,T54 

,'CO 

— 

', 

•  ,16 

,T28,' INCREMENT  = 

«, 

13, 

,T28, 

•INCREMENT  = 

«,I3 

rT54, 

,15X, 
MR!  I) 

•B' ,9X,'MAX. 
,I=i,NS) 

RADIUS' 

) 

WZI=WC/ITEMPR 

DI=PIC*(WZI**0.322) 

WRAT*WZI/(WW*2 .0*P I *1 OOO.O) 

WRAT2=WRAT*WRAT 

WRATM=WRAT2-1.0 

BATH=WRATM/( { WRATM*WRATM)+< DI*D I*WRAT2  )  ) 
C    SELECT    U(R)D(R). 

DO    920    J=1,NS 

IF  (ITEMPR  .GT.MR(J))  GO  TO  920 

URDR=1.0E-09*( A(J)*ITEMPR+B(J) ) 

BUBBLE=URDR/( ITEMPR* ITEMPR*! .OE-12) 

GO  TO  927 
920  CONTINUE 
C  IF  FALLS  THROUGH  HERE,  WE  HAVE  GONE  TOO  FAR. 

925  WRITE(6,926) 

926  FORMAT(»  ERROR  ON  SEGMENT  RADIUS.') 
STOP 

C  WRITE  RESULT  FOR  EACH  I. 

927  IF  (K.NE.l)  GO  TO  940 
WRITE(6,930)  L , ITE MPR , WZI ,URDR 

9  30  FORMAT (2X, 16 , 6X, 16 , 6X, E16. 8, 2X , E16.8I 
940  ITEMPR=ITEMPR+IDR 

SUM=SUM+(BUBBLE*BATH) 
1000  CONTINUE 
C  NOW  CALCULATE  SOUND  VELOCITY. 

CNEW=1.0/(  (1.0/CZ)+U3.0*CZ*SUM)/(8,0*PI*PI*WW*WW* 
11.0E06) ) ) 
CC(K)=CNEW-CZ 

WRITE ( 6,1010)  ITEMPW,SUM,CNEW,CC(K) 
1010  FORMAT('  W  =  ,,I6,4X,«SUM  =  • , E16. 8,2X, »C  =  «,E16.8, 
12X,'DIFF  =  • ,E16.8  ) 
C  PREPARE  FOR  NEXT  ROUND  OF  W. 
ITEMPW=ITEMPW+IDW 
2000  CONTINUE 
C  ALL  DONE.. NOW  PLOT  A    AND  C. 
WRITE (6, 21 00) 
2100  FORMAT(»l  ^FSULT:   W  VS.  (C-CO)») 

WRITE< 6,2200 )( W(I ) ,CC( I ) ,I=1,IW) 
22  00  F0RMAT(10X,E20.10, 10X, E20.10) 

WRITE(6,2300) 
2300  FORMAT ( • 1« ) 

CALL  PLOTP  (W,CC,IW,0) 

STOP 

END 
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